Simulating and Trading Endogenous Feedback in the Stock Market
Modelling Stock Price Dynamics, Implied Volatility and Black-Scholes Derivatives
Gary He
1 Endogenous Feedback Model
We set out to simulate an endogenous market feedback model as laid out in Danielsson, Shin and Zigrand’s ‘Endogenous Extreme Events and the Dual Role of Prices’ (2012). In light of the recent craze in retail option investing, which has driven market makers, hedge funds and short sellers to the brink of collapse whilst creating millionaires amongst amateur traders in the space of a few months, we are motivated to understand how the hedging of option positions can create feedback loops in prices.
A simple toy model for the stock price dynamics with feedback generated by market makers is given as a variation of the standard geometric Brownian motion in Danielsson et al’s paper:
\[ dS_t = {\mu S_t dt + \sigma S_t dW_t - \psi_t d\Delta_t} \]
\[ \mu \geq 0 \text{ and } \sigma \geq 0 \text{ and }\psi_t := \psi S_t \theta \text{ or } \psi_t \theta \] The intuition behind this model is simple. We have a impact term, psi, which permanently pushes prices up or down as delta, the change in option price for change in the underlying, shifts. Assuming that market makers are delta hedging options they have sold into the market, they will be short delta and short gamma. As the stock price goes up, market makers will have to buy more of the underlying to hedge, thereby increasing the stock price further. The same applies for the reverse, giving us a simple, one-dimensional feedback loop model for the stock price.
Using some stochastic calculus and Girsanov’s change of measure, we can rearrange this expression into a standard geometric Brownian motion with drift, diffusion, and no extraneous terms.
We apply Ito’s lemma to delta, assuming that delta is a function of time and the underlying:
\[ d\Delta_t = \frac{d\Delta_t}{dt}dt + \frac{d\Delta_t}{dS_t}dS_t+\frac{1}{2}\frac{d^2\Delta_t}{dS_t^2}dS_t^2 \] \[ \text{ where } dS_t^2 = \sigma^2 S_t^2 dt \text{ and } \frac{d^2P}{dS_t^2}=\frac{d\Delta_t}{dt}:=\Gamma_t \] \[ d\Delta_t = \frac{d\Delta_t}{dt}dt + \frac{d\Delta_t}{dS_t}dS_t+\frac{1}{2}\frac{d^2\Delta_t}{dS_t^2}\sigma^2S_t^2dt \] As we are able to express higher order differentials of the stock price in terms of Gamma, Charm and Speed, we can simplify the expression further:
\[ \frac{d\Delta_t}{dt} = -\phi (d_1) \frac{2r \tau - d_2 \sigma \sqrt{\tau}}{2 \tau \sigma \sqrt{\tau}} = \frac{-\Gamma_t S_t (2 r \tau - d_2 \sigma \sqrt{\tau})}{2 \tau} \]
\[ \frac{d^2\Delta_t}{dS_t^2} = -\frac{\Gamma_t}{S_t} (\frac{d_1}{\sigma \sqrt{\tau}}+1) \]
\[ \frac{1}{2} \frac{d^2\Delta_t}{dS_t^2} \sigma^2 S_t^2 = -\frac{\Gamma_t \sigma S_t (d_1 + \sigma \sqrt{\tau})}{2 \sqrt{\tau}} \]
\[ d\Delta_t = \frac{d\Delta_t}{dt}dt + \frac{d\Delta_t}{dS_t}dS_t+\frac{1}{2}\frac{d^2\Delta_t}{dS_t^2}\sigma^2S_t^2dt = -\Gamma_t S_t (r + \sigma^2)dt + \Gamma_t dS_t \]
We now bring in Girsanov’s change of measure. Under the risk neutral measure where the stock price, rather than the risk-free money market instrument is used as a numeraire, we can show that the real-world drift term corresponds to the following:
\[ \mu = (r + \sigma^2) \] Combining these ‘instantaneous adjustments’, we obtain the following basic expression for the change in delta:
\[ d\Delta_t = -\Gamma_t \mu S_t dt + \Gamma_t dS_t \]
By substituting in our new change in delta expression to the original formula, we can arrive at a standard geometric Brownian motion with drift and diffusion as promised: \[ dS_t = \mu S_t dt + \sigma S_t dW_t - \psi_t[-\Gamma_t \mu S_t dt + \Gamma dS_t] \]
\[ dS_t = \mu S_t dt + \frac{\sigma S_t}{1 + \psi_t \Gamma_t} dW_t \] This form of the stock price dynamics will form the basis for our market simulations, as we can use Monte Carlo simulations, discretisation procedures, and our understanding of options to analyse how stocks behave under this hypothetical scenario.
Firstly, we will create some basic functions for programming the option ‘Greeks’, which will be required to calculate the actual stock price. These functions will also be important for tracking the behaviour of options and implied volatility over time, just as practitioners and traders do in reality.
2 Option Dynamics
2.1 First-Order Greeks
We begin with the simple first-order Greeks, which are simply the first derivatives of the option parameters. We also calculate the price of puts/calls and implied volatility, as we particularly would like to understand how volatility can be affected in these market situations.
2.1.1 Function Parameters
Here we set out some basic parameters that are common across our option formulas. We also load up necessary libaries for data visualisation:
#pc put/call indicator call=1, put=-1
#S stock price at 0
#K strike
#price option premium
#d dividend yield
#r risk-free rate
#t time to maturity
#start starting value for Newton method estimation for implied volatility, optional, by default=0.2
library(ggplot2)
library(plotly)2.1.2 Price
We derive the price from the famous Black-Scholes model, which also forms the basis for the Greeks:
\[ V = \pm(Se^{-dt} \Phi(d_1) - Ke^{-rt} \Phi(d_2)) \] \[ \text{where } d_1 = \frac{ln\frac{S}{K}+(r-d+\frac{\sigma^2}{2})}{\sigma \sqrt{\tau}} , d_2 = d_1 - \sigma \sqrt{\tau} \text{ and } \tau=T-t \]
BSprice<-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSprice = pc * exp(-d * t) * S *
pnorm(pc * d1) - pc * k * exp(-r * t) * pnorm(pc * d2)
return(BSprice)
}2.1.3 Delta
Delta is defined as the change in option value for change in the underlying, and is fundamental to linear hedging.
\[ \Delta = \frac{dV}{dS} = e^{-dt} \Phi(d_1) \]
BSdelta<-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
if (pc == 1) {BSdelta = exp(-d * t) * pnorm(d1)} else
{BSdelta = exp(-d * t) * (pnorm(d1) - 1)}
return(BSdelta)
}2.1.4 Gamma
The second derivative of the option value with respect to the underlying, Gamma, measures convexity in option value.
\[ \Gamma = \frac{d\Delta}{dS} = e^{-dt} \Phi(d_1) \]
BSgamma<-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
BSgamma = exp(-d * t) * exp((-d1 ^ 2) / 2) / (sqrt(2 * pi) * S * vol * sqrt(t))
return(BSgamma)}2.1.5 Theta
Theta is the derivative with respect to time-to-maturity. Unsurprisingly, it is usually negative so as to satisfy no arbitrage conditions.
\[ \theta = \frac{dV}{d\tau} = -re^{-d\tau}K\Phi(d_2)-\frac{S\sigma}{2\sqrt{\tau}}\phi(d_1) \] \[ \text{where } \phi(d_1) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}d_1^2} \]
BStheta<-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BStheta = -exp(-d * t) * exp((-d1 ^ 2) / 2) * S * vol /
(sqrt(2 * pi) * 2 * sqrt(t)) + pc * d * S * exp(-d * t) * pnorm(pc * d1) - pc * r * k * exp(-r * t) * pnorm(pc * d2)
return(BStheta)
}2.1.6 Vega
Vega is the derivative with respect to volatility. Contrary to theta, it is strictly positive, as volatility tends to increase likelihood of ending ITM.
\[ \nu = \frac{dV}{d\sigma} = S\phi(d_1)\sqrt{\tau} \]
BSvega<-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
BSvega = exp(-d * t) * S * sqrt(t) * exp((-d1 ^ 2) / 2) / (sqrt(2 * pi))
return(BSvega)
}2.1.7 Rho
Rho represents the derivative with respect to the risk-free interest rate.
\[ \rho = \frac{dV}{dr}\pm KTe^{-rT}\Phi(\pm d_2) \]
BSrho<-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSrho = pc * k * t * exp(-r * t) * pnorm(pc * d2)
return(BSrho)
}2.1.8 Implied volatility
Implied volatility is a vital tool for us to price options within the market model. It can be seen as the solution to the following problem:
\[ \text{solve } V(S,K,r,\tau,\sigma) = V_{BS}(S,K,r,\tau,\sigma_{IV}) \text{ for } \sigma \]
To do this, we can either use unit roots as programmed in R:
BSvol1 <- function(pc, S, k, d, r, t, price)
{
root_fun = function(vol) {
BSprice(pc,S,k,vol,d,r,t) - price
}
uniroot(root_fun,c(0,10))$root
}Or we may use Newton’s gradient method:
BSvol<-function(pc, S, k, price, d, r, t, start = 0.2)
{
voli = start
pricei = BSprice(pc, S, k, voli, d, r, t)
vegai = BSvega(pc, S, k, voli, d, r, t)
while(abs(price - pricei) > 0.000001)
{
voli<-voli + (price - pricei) / vegai
pricei<-BSprice(pc, S, k, voli, d, r, t)
vegai<-BSvega(pc, S, k, voli, d, r, t)
}
BSvol = voli
return(BSvol)
}2.2 Higher-Order Option Greeks
We also include higher-order Greeks, which tend to play a large role in FX markets, though are less commonly given attention in equity trading.
2.2.1 Charm
\[ Charm = -\frac{d^2V}{dSd\tau} = -\phi(d_1)\frac{2r\sqrt{\tau}-\sigma d_2}{2\sigma\tau} \]
BScharm <-function(pc, S, k, vol, d, r, t)
{
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
if (pc == 1) {
BScharm = -1*(-(1/(sqrt(2*pi)))*exp(-d1^2/2)*(2*r*sqrt(t)-vol*d2)/(2*vol*t))} else
{BScharm = -(1/(sqrt(2*pi)))*exp(-d1^2/2)*(2*r*sqrt(t)-vol*d2)/(2*vol*t)}
return(BScharm)
}2.2.2 Vanna
\[ Vanna = \frac{d^2V}{d\sigma dS} = \sqrt{\tau} \phi(d_1) [1-d_1] \]
Version 1:
BSvanna<-function(pc, S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSvanna2 = -BSgamma(pc, S, k, vol, d, r, t) * (sqrt(t) * S) * d2
return(BSvanna2)
}Version 2:
BSvanna2<-function(pc, S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSvanna3 = (BSvega(pc, S, k, vol, d, r, t) / S) * (1 - d1 / (vol * sqrt(t)))
return(BSvanna3)
}2.2.3 VegaVanna
\[ VegaVanna = \frac{d^3V}{d\sigma^2 dS} \]
BSvegavanna<-function(pc, S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSvegavanna2 = BSvanna2(pc,S,k,vol,d,r,t) * (1 / vol) * (d1 * d2 - d1 / d2 - 1)
return(BSvegavanna2)
}2.2.4 Veta
\[ Veta = \frac{d^2V}{d\sigma d\tau} \]
BSveta<-function(S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSveta = -S * exp(-d * t) * exp((-d1 ^ 2) / 2) / (sqrt(2 * pi)) * sqrt(t) * (d + (((r-d) * d1) / (vol * sqrt(t))) - ((1 + d1 * d2) / 2*t))
return(BSveta)
}2.2.5 Volga
\[ Volga = \frac{d^2V}{d\sigma^2 } \]
BSvolga<-function(S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSvolga = exp(-d * t) * sqrt(t) * exp((-d1 ^ 2) / 2) / (sqrt(2 * pi)) * ((d1 * d2) / vol)
return(BSvolga)
}2.2.6 Speed
\[ Speed = \frac{d^3V}{dS^3} \]
BSspeed<-function(pc, S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSspeed = (-BSgamma(pc, S, k, vol, d, r, t) / S) * ((d1 / (vol * sqrt(t))) + 1)
return(BSspeed)
}2.2.7 Zomma
\[ Zomma = \frac{d^3V}{d\sigma dS^2} \]
BSzomma<-function(pc, S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSzomma = BSgamma(pc, S, k, vol, d, r, t) * ((d1 * d2 - 1) / vol)
return(BSzomma)
}2.2.8 Ultima
\[ Ultima = \frac{d^2S}{d\sigma^3 } \]
BSultima<-function(pc, S, k, vol, d, r, t){
d1 = (log(S / k) + t * (r - d + (vol ^ 2) / 2)) / (vol * sqrt(t))
d2 = d1 - vol * sqrt(t)
BSultima = (-BSvega(pc, S, k, vol, d, r, t) / vol^2) * (((d1 * d2) * (1 - d1 * d2)) + d1^2 + d2^2)
return(BSultima)
}3 Simulation - Model 1
We now are equipped with the tools and functions required to begin simulating stock price paths. To do this, we can use the Milstein discretisation scheme, which has greater accuracy compared to other methods, such as the Euler method, in approximating continuous functions. As the change in a standard Brownian motion is normally distributed, we can also use the Box-Muller transform to generate random i.i.d numbers. We choose the following parameters for our model:
#pc = 1 #market makers have sold call options into the market
#mu = 0 #drift is 0 under the risk neutral model
#r = 0 #discount rate is 0
#dt = 1/12000 #we split a month into 1000 step
#nsteps = 1000 #number of discretised periods
#stock = 100 or 120 #initial stock price
#theta = 1 #number of options being hedged
#psi = -0.1 #permanent impact on prices is negative due to market makers short net positions
#K = 100 or 120 #market makers begin by hedging ATM
#ttm = 1/12 #maturity of market makers options is a month
#strike = #range of strikes will be used for volatility skewsWe can now begin generating stock price paths, which predictably behave like a geometric Brownian motion:
getpricepath = function(pc,mu,vol,dt,nsteps,stock,theta,psi,K,ttm) {
U1<-runif(nsteps)
U2<-runif(nsteps)
N_0_1<-sqrt(-2*log(U1))*cos(2*pi*U2)
W<-c()
W<-N_0_1*sqrt(dt)
stock_price <<- c()
stock_price[1] <<- stock
for (i in 2:length(W))
{
stock_price[i]<<-stock_price[i-1] + mu * stock_price[i-1] * dt + ((vol) / (1 + theta * psi * stock_price[i-1] * BSgamma(pc,stock_price[i-1],K,vol,0,0,ttm-(i-2)*dt))) * stock_price[i-1] * W[i-1]+ 0.5 * (W[i-1]^2 - dt) * stock_price[i-1] * ((vol)/(1 + theta * psi * stock_price[i-1] * BSgamma(pc,stock_price[i-1],K,vol,0,0,ttm-(i-2)*dt)))^2
}
plot(stock_price,type="l",xlab="Steps",ylab="Price")
}
getpricepath(1,0,0.3,1/12000,1000,100,-0.1,1,100,1/12)
Alternatively, we can use the built in cumsum and rnorm functions to reduce the computational load of generating stock price paths, though this fails to capture the changing gamma of market makers over time. This procedure could be useful however for a discrete hedging model. For illustrative purposes however, we assume that market makers will hedge constantly.
We wish to capture moderate effects of market feedback as per Danielsson et al’s paper, such that the Brownian coefficient \[ \frac{\sigma S_t}{1 + \psi_t \Gamma} \] is equal to 2. We find that to do this, we can set the the product of psi and theta to -0.1 at initiation. We expect implied volatility to be significantly amplified at this level of specified feedback, as the impact of Brownian shocks is effectively doubled.
To see how the stock price behaves, we can generate and plot simulated price paths:
getpricepath1 = function(pc,mu,vol,dt,nsteps,stock,theta,psi,K,ttm) {
x<-(mu-vol^2/2)*dt + vol*sqrt(dt)*rnorm(nsteps,0,1)*((vol) / (1 + theta * psi * BSgamma(pc,stock,K,vol,0,0,1/12)))
price = stock * exp(cumsum(x))
plot(price, type="l", xlab="Steps", ylab="Price")
}
getpricepath1(1,0,0.3,1/12000,1000,100,-0.1,1,100,1/12)
For functional purposes however, we simply wish to obtain the final price of a simulated path. This will allow us to simulate reasonable option prices and implied volatilities:
getprice = function(pc,mu,vol,dt,nsteps,stock,theta,psi,K,ttm) {
U1<-runif(nsteps)
U2<-runif(nsteps)
N_0_1<-sqrt(-2*log(U1))*cos(2*pi*U2)
W<-c()
W<-N_0_1*sqrt(dt)
stock_price <<- c()
stock_price[1] <<- stock
for (i in 2:length(W))
{
stock_price[i]<<-stock_price[i-1] + mu * stock_price[i-1] * dt + ((vol) / (1 + theta * psi * stock_price[i-1] * BSgamma(pc,stock_price[i-1],K,vol,0,0,ttm-(i-2)*dt))) * stock_price[i-1] * W[i-1]+ 0.5 * (W[i-1]^2 - dt) * stock_price[i-1] * ((vol)/(1 + theta * psi * stock_price[i-1] * BSgamma(pc,stock_price[i-1],K,vol,0,0,ttm-(i-2)*dt)))^2
}
return(stock_price[nsteps])
}
getprice(1,0,0.3,1/12000,1000,100,-0.1,1,100,1/12)## [1] 111.46483.1 Implied Volatility
More importantly however, we want to understand how implied volatility is impacted using this model. We build a simple function to return the implied volatility of a given option, with specified time to maturity and strike, within our market model:
mcvol = function(nreps,pc,mu,vol,dt,nsims,stock,theta,psi,K,ttm,strike) {
c=c()
for(i in 1:nreps) {
St = getprice(pc,mu,vol,dt,nsims,stock,theta,psi,K,ttm)
c[i] = max(St-strike,0)
}
c = na.omit(c)
return(BSvol(pc,stock,K,mean(c),0,0,nsims*dt))
}
mcvol(1000,1,0,0.3,1/12000,500,100,-0.1,1,100,1/12,100)## [1] 0.56458193.1.1 Volatility Skews Over Time
We would like to see whether our market model creates a ‘realistic’ volatility smirk as observed with real-world equity markets. As described by the so-called ‘leverage effect’, we would expect there to be a positive skew of volatility as strikes decrease. Furthermore, we should expect volatility to flatten over time, in-line with the limitations of the Black-Scholes model. We can first create a 2D superposed line graph of the volatility smirk at each time period. We do this first, assuming that the starting price is 120, the simulation runs for a month, that market makers are hedging an initially ATM option.
ivplot = function(nreps,range1,range2,interval,pc,mu,vol,dt,nsims,stock,theta,psi,K,ttm) {
s = c()
s = seq(range1,range2,interval)
v = c()
v2 = c()
v3 = c()
v4 = c()
v5 = c()
v6 = c()
v7 = c()
v8 = c()
v9 = c()
v10 = c()
for(i in 1:length(s)) {
v[i] = mcvol(nreps,pc,mu,vol,dt,nsims,stock,theta,psi,K,ttm,s[i])
v2[i] = mcvol(nreps,pc,mu,vol,dt,2*nsims,stock,theta,psi,K,ttm,s[i])
v3[i] = mcvol(nreps,pc,mu,vol,dt,3*nsims,stock,theta,psi,K,ttm,s[i])
v4[i] = mcvol(nreps,pc,mu,vol,dt,4*nsims,stock,theta,psi,K,ttm,s[i])
v5[i] = mcvol(nreps,pc,mu,vol,dt,5*nsims,stock,theta,psi,K,ttm,s[i])
v6[i] = mcvol(nreps,pc,mu,vol,dt,6*nsims,stock,theta,psi,K,ttm,s[i])
v7[i] = mcvol(nreps,pc,mu,vol,dt,7*nsims,stock,theta,psi,K,ttm,s[i])
v8[i] = mcvol(nreps,pc,mu,vol,dt,8*nsims,stock,theta,psi,K,ttm,s[i])
v9[i] = mcvol(nreps,pc,mu,vol,dt,9*nsims,stock,theta,psi,K,ttm,s[i])
v10[i] = mcvol(nreps,pc,mu,vol,dt,10*nsims,stock,theta,psi,K,ttm,s[i])
}
voll<<-cbind(v,v2,v3,v4,v5,v6,v7,v8,v9,v10)
}
ivplot(100,90,130,5,1,0,0.3,1/12000,100,120,-0.1,1,120,1/12)
voll1 <- data.frame(voll)
voll120 <- voll
strikes = seq(90,130,5)
vol = voll1[,1]
ggplot(voll1,mapping=aes(x=strikes,y=vol)) + geom_line(aes(y=v),color="black") + geom_line(aes(y=v2),color="red") + geom_line(aes(y=v3),color="blue") + geom_line(aes(y=v4),color="brown") + geom_line(aes(y=v5),color="skyblue") + geom_line(aes(y=v6),color="slateblue") + geom_line(aes(y=v7),color="orange") + geom_line(aes(y=v8),color="tomato") + geom_line(aes(y=v9),color="olivedrab4") + geom_line(aes(y=v10),color="mistyrose2") + theme(panel.background = element_blank(), panel.border = element_rect(color = "black",fill=NA), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + labs(title="Implied Volatility 2D Surface") 
We do the same with the assumption that the starting price is instead 100, to see if we can gain any insights into the trader regime in our market model, as proposed by Derman.
ivplot(100,90,130,5,1,0,0.3,1/12000,100,100,-0.1,1,100,1/12)
voll2 <- data.frame(voll)
voll100 <- voll
strikes = seq(90,130,5)
vol = voll2[,1]
ggplot(voll2,mapping=aes(x=strikes,y=vol)) + geom_line(aes(y=v),color="black") + geom_line(aes(y=v2),color="red") + geom_line(aes(y=v3),color="blue") + geom_line(aes(y=v4),color="brown") + geom_line(aes(y=v5),color="skyblue") + geom_line(aes(y=v6),color="slateblue") + geom_line(aes(y=v7),color="orange") + geom_line(aes(y=v8),color="tomato") + geom_line(aes(y=v9),color="olivedrab4") + geom_line(aes(y=v10),color="mistyrose2") + theme(panel.background = element_blank(), panel.border = element_rect(color = "black",fill=NA), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) + labs(title="Implied Volatility 2D Surface") 
3.1.2 3D Volatility Surface
We can refine our findings further by creating a 3D implied volatility surface, which also lets us view the term structure of implied volatility in a dynamic fashion. Superimposing multiple surfaces on a single graph also allows us to assess which trading regime we are in, by combining the generated 2D skews for different values of the underlying stock:
strikes = seq(90,130,5)
ttm = seq(0,30,3)
fig<- plot_ly(x=strikes, y=ttm, showscale = FALSE)
fig<- fig %>% add_surface(colorscale = list(seq(0,1,length.out=7),c("#18389B","#02B1E8","#02B283","#FEBB08","#F93616","#F7103C","#CC0B8C")),z=~voll100)
fig<- fig %>% add_surface(colorscale = list(seq(0,1,length.out=7),c("#B2182B","#D6604D","#F4A582","#FDDBC7","#D1E5F0","#92C5DE","#4393C3")),z=~voll120)
fig<- fig %>% layout(
title = "Implied Volatility Surface",
scene = list(
xaxis = list(title = "Strike"),
yaxis = list(title = "Time to Maturity (Days)"),
zaxis = list(title = "Implied volatility (%)",tickformat="%"),
aspectratio = list(x = 1, y = 1, z = 1),
camera = list(eye =list(x=1.25,y=1.25,z=0.5))
))
figWe see evidence that we may be in a ‘sticky delta’ regime, as volatility changes as a function of the underlying and is not constant with the strike. Note that due to limitations in computational power, we only run 100 simulations at each strike, which may lead to an uneven surface. With a higher number of simulations, the surfaces are more consistent.
4 Simulation - Model 2
We now consider a more refined model in which the drift is determined by the expected delta hedging by market makers over time, given that Riemann integrals are equivalent to expectations:
\[ dS_t = - \psi M_t^{\Delta}dt + \sigma S_t dW_t \text{ } ;M_t^{\Delta}dt:=\int_0^t e^{\varphi (t-u)} d\Delta_u \] We will simulate stock prices as before, but will adjust parameters slightly to remain consistent with the previous model:
#pc = 1 #market makers have sold call options into the market
#mu = 0 #drift is 0 under the risk neutral model
#r = 0 #discount rate is 0
#dt = 1/12000 #we split a month into 1000 step
#nsteps = 1000 #number of discretised periods
#stock = 100 or 120 #initial stock price
#theta = -0.3 #number of options being hedged
#psi = -0.4 #permanent impact on prices is negative due to market makers short net positions
#K = 100 or 120 #market makers begin by hedging ATM
#ttm = 1/12 #maturity of market makers options is a month
#strike = #range of strikes will be used for volatility skews4.1 Tracking Option Dynamics
We can also track over time the first and higher order Greeks of options where the stock is the underlying asset, which may provide us with insight into the price dynamics and behaviour of stocks in our simulated model:
getgreeks2 = function(pc,S,k,vol,nsteps,dt,psi,theta,a,div,r) {
N_0_1<-rnorm(nsteps)
W<-c()
W<-N_0_1*sqrt(dt)
stock_price<- term1 <- d <- g <- spd <- van <- chm <- vet <- vog <- c()
stock_price[1] = 0
stock_price[2] = S
d[1] = 0
d[2] = BSdelta(pc,S,k,vol,div,r,dt*nsteps)
term1[1] = 0
term1[2] = psi*theta*S*exp(-dt)*(d[2]-d[1])*dt
g[1] = 0
g[2] = BSgamma(pc,S,k,vol,div,r,dt*nsteps)
for (i in 3:length(W))
{
term1[i] = exp(-dt*a)*(d[i-1]-d[i-2])*dt
stock_price[i]<-stock_price[i-1]+psi*theta*stock_price[i-1]*sum(term1) + vol*stock_price[i-1]*W[i-1]+0.5*(W[i-1]^2-dt)*stock_price[i-1]*vol^2
d[i] = BSdelta(pc,stock_price[i],k,vol,div,r,(nsteps-i)*dt)
g[i] = BSgamma(pc,stock_price[i],k,vol,div,r,(nsteps-i)*dt)
spd[i] = BSspeed(pc,stock_price[i],k,vol,div,r,(nsteps-i)*dt)
van[i] = BSvanna(pc,stock_price[i],k,vol,div,r,(nsteps-i)*dt)
chm[i] = BScharm(pc,stock_price[i],k,vol,div,r,(nsteps-i)*dt)
vet[i] = BSveta(stock_price[i],k,vol,div,r,(nsteps-i*dt))
vog[i] = BSvolga(stock_price[i],k,vol,div,r,(nsteps-i*dt))
}
g <- g[-(1)]
stock_price <- stock_price[-(1)]
d <- d[-(1)]
par(mfrow=c(2,2))
plot(stock_price,type="l")
plot(g,type="l",col="blue")
plot(d,type="l",col="red")
plot(spd,type="l",col="orange")
plot(van,type="l",col="purple")
plot(chm,type="l",col="green")
plot(vet,type="l",col="yellow")
plot(vog,type="l",col="pink")
}
getgreeks2(1,100,90,0.3,1000,1/12000,-0.3,-0.4,1,0,0)
We can gather that in such a model, as delta increases, feedback loops cause realised volatility to increase, such that it likely exceeds implied volatility. Conversely, as delta decreases, volatility is likely to decrease, approaching implied volatility.
4.2 Algorithmic Trading Strategy 1: Gamma Scalping
One way of trading this would simply to Gamma scalp - a delta neutral strategy whereby going long an option and adjusting a short position in the underlying over time, should allow us to earn realised volatility whilst shorting implied. Arguably, in such a market model, there is ample opportunity for realised volatility to exceed implied due to the feedback effect. We construct a simple function to carry out this trade in real-time, rather than depending on a stored price path, to see whether this strategy is profitable. We plot the resulting distribution:
trade_1 = function(pc,S,k,vol,nsteps,dt,psi,theta,phi,div,r) {
N_0_1 <- rnorm(nsteps)
W <- c()
W <- N_0_1 * sqrt(dt)
stock_price <- term1 <- d <- price <- gain <- c()
stock_price[1] <- term1[1] <- d[1] <- price[1] <-gain[1:2]<- 0
stock_price[2] = S
d[2] = BSdelta(pc,S,k,vol,div,r,dt * nsteps)
term1[2] = psi * theta * S * exp(-phi * dt) * (d[2]-d[1])*dt
price[2] = BSprice(pc,S,k,vol,div,r,dt*nsteps)
pl <- rep(0, length(stock_price))
pl[2] <- price[2]+d[2]*stock_price[2]
for (i in 3:length(W))
{
term1[i] <- exp(-phi * dt)*(d[i-1]-d[i-2])*dt
stock_price[i] <- stock_price[i-1] + psi * theta * stock_price[i-1] * sum(term1) + vol * stock_price[i-1] * W[i-1] + 0.5 * (W[i-1]^2-dt) * stock_price[i-1] * vol^2
d[i] <- BSdelta(pc,stock_price[i],k,vol,div,r,(nsteps - i) * dt)
price[i] <- BSprice(pc,stock_price[i],k,vol,div,r,(nsteps - i) * dt)
pl[i] <- d[i] * stock_price[i]
gain[i] <- pl[i] - pl[i-1]
}
return(((pl[2] + sum(gain)) / pl[2]) - 1) }
trade_1(1,100,100,0.3,1000,1/12000,-0.3,-0.4,1,0,0)## [1] 1.112701trade_1_dist = c()
for(i in 1:1000) {
trade_1_dist[i] = 100*trade_1(1,100,100,0.3,1000,1/12000,-0.3,-0.4,1,0,0)
}
hist(trade_1_dist, col="red", border="darkred",xlab="Returns in %", main="Trade 1 PnL", prob = TRUE,breaks=50)
lines(density(trade_1_dist), lwd = 2, col = "blue")
summary(trade_1_dist)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -100.00 82.13 93.76 53.89 106.18 158.88It is clear that Gamma scalping can be profitable, likely when Gamma is high and realised volatility is also high. However, we see that this strategy carries significant tail risk - we may completely lose our investment due to the failure of feedback loops to generate significant realised volatility.
4.3 Algorithmic Trading Strategy 2: Feedback Capture and Gamma Scalping
Instead, we can combine our understanding of the delta feedback loop with our knowledge of how Gamma scalping can allow us to go long realised and short implied volatility. We construct a trading algorithm that will go long realised as delta increases, and short realised as delta decreases, attempting to capture the feedback loop. We also implement a buy discipline that will go short and hold as the option moves far ITM such that delta is 1. At this point, the volatility will have peaked. We plot the resulting P&L distribution:
trade_2 = function(pc,S,k,vol,nsteps,dt,psi,theta,phi,div,r) {
N_0_1 <- rnorm(nsteps)
W <- c()
W <- N_0_1 * sqrt(dt)
stock_price <- term1 <- d <- price <- c()
stock_price[1] <- term1[1] <- d[1] <- price[1] <- 0
stock_price[2] <- S
d[2] <- BSdelta(pc,S,k,vol,div,r,dt * nsteps)
term1[2] = psi * theta * S * exp(-phi * dt) * (d[2] - d[1]) * dt
price[2] = BSprice(pc,S,k,vol,div,r,dt*nsteps)
pl <- rep(0, length(stock_price))
longpl <- rep(0, length(stock_price))
shortpl <- rep(0, length(stock_price))
for (i in 3:length(W))
{
term1[i] <- exp(-phi * dt) * (d[i-1] - d[i-2]) * dt
stock_price[i] <- stock_price[i-1] + psi * theta * stock_price[i-1] * sum(term1) + vol * stock_price[i-1] * W[i-1] + 0.5 * (W[i-1]^2 - dt) * stock_price[i-1] * vol^2
d[i] <- BSdelta(pc,stock_price[i],k,vol,div,r,(nsteps - i) * dt)
price[i] <- BSprice(pc,stock_price[i],k,vol,div,r,(nsteps - i) * dt)
longpl[i] <- -price[i] + (d[i] * stock_price[i])
shortpl[i] <- -(d[i] * stock_price[i]) + price[i]
}
sell = FALSE
buy = TRUE
trade <- rep(0,length(stock_price))
for(i in 2:length(stock_price)) {
if(d[i] < d[i-1] && buy && !isTRUE(all.equal(d[i],1))) {
sell = TRUE
buy = FALSE
trade[i] = 2
} else if(d[i] > d[i-1] && sell) {
trade[i] = 1
sell = FALSE
buy = TRUE
} else if(d[i] < d[i-1] && sell) {
trade[i] = 0
sell = TRUE
buy = FALSE
} else if(d[i] > d[i-1] && buy) {
trade[i] = 0
buy = TRUE
sell = FALSE
} else if(isTRUE(all.equal(d[i],1)) && buy || sell) {
trade[i] = 2
buy = FALSE
sell = FALSE
} else if(isTRUE(all.equal(d[i],1)) && !buy && !sell) {
trade[i] = 0
}
}
for(i in 2:length(stock_price)) {
if(trade[i] == 1) {
pl[i] = longpl[i]
} else if(trade[i] == 2) {
pl[i] = shortpl[i]
} else if(trade[i] == 0) {
pl[i] = 0
}
}
row = which(trade!=0)
if(trade[length(stock_price)]==0) {if(trade[row[length(row)]] ==1) {pl[length(stock_price)] = shortpl[length(stock_price)]} else if(trade[row[length(row)]] ==2) {pl[length(stock_price)] = longpl[length(stock_price)]}}
return((sum(pl)+abs(pl[row[1]]))/(abs(pl[row[1]])))
}
trade_2(1,100,100,0.3,1000,1/12000,-0.3,-0.4,1,0,0)## [1] 2.087458trade_2_dist = c()
for(i in 1:1000) {
trade_2_dist[i] = 100*trade_2(1,100,100,0.3,1000,1/12000,-0.3,-0.4,1,0,0)
}
hist(trade_2_dist, col="red", border="darkred",xlab="Returns in %", main="Trade 2 PnL", prob = TRUE, breaks = 100)
lines(density(trade_2_dist), lwd = 2, col = "blue")
summary(trade_2_dist)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -197.36 62.38 130.19 109.13 164.28 305.09We can see that this strategy is highly profitable, with significantly lower tail risk than purely Gamma scalping which still contains directional risk. Whilst the strategy is not flawless, it accurately captures the behaviour of volatility and delta hedgers and their impact on stock prices.